Viscoelasticity of Dilute Solutions of Semiflexible Polymers 
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We show using Brownian dynamics simulations and theory how the shear relaxation modulus G(t) 
of dilute solutions of relatively stiff semiflexible polymers differs qualitatively from that of rigid rods. 
For chains shorter than their persistence length, G(t) exhibits three time regimes: At very early 
times, when the longitudinal deformation is affine, Git) ~ t~ 3 ^ 4 . Over a broad intermediate regime, 
during which the chain length relaxes, G{i) ~ t~ 5 ^ 4 . At long times, G(t) mimics that of rigid rods. 
A model of the polymer as an effectively extensible rod with a frequency dependent elastic modulus 
B(ui) ~ (ilj) 3//4 quantitatively describes G(t) throughout the first two regimes. 



Many important biopolymers are wormlike chains with 
persistence lengths L p comparable to or larger than their 
contour length L. Examples are a-helical proteins, short 
DNA, collagen fibrils, rod-like viruses, and protein fila- 
ments such as F-actin. The cytoskeleton of a cell is pri- 
marily a network of such polymers, and plays a critical 
role in controlling the mechanical rigidity, motility, and 
adhesion of living cells; understanding the viscoelastic 
behavior semiflexible polymers in solution is thus a criti- 
cal problem in biophysics. Whereas the linear viscoelas- 
tic behavior of dilute solutions of flexible (Gaussian) and 
rod- like polymer molecules is well understood JjJ, there 
is thus far no qualitatively correct description of the vis- 
coelasticity of dilute solutions of semilexible polymers 
over the whole range of frequency and time scales. Bridg- 
ing the theoretical gap between the flexible and rigid rod 
limits is thus also an important open problem in polymer 
physics. Here, we present both results from Brownian dy- 
namics simulations of relatively stiff semiflexible chains, 
with L p > L, and a simple theory that accurately de- 
scribes their linear viscoelastic response over a very wide 
range of time scales. Both theory and simulation yield 
a relaxation modulus G(t) ~ i~ 5 / 4 over a wide range of 
intermediate times, after an initial decay of Git) oc i -3 / 4 
at very early times, and before an exponential decay of 
Git) at long times (like that of rigid rods) due to diffusive 
tumbling of the chain orientation. 

A single wormlike chain may be described by a curve 
r(s), with a tangent vector u(s) = dr(s)/ds, where 
s is contour distance along the chain. Inextensibility 
requires that \dr/ds\ = 1. The bending energy of a 
chain with rigidity k or persistence length L p = n/kT 
is U = ^njds | | 2 . The Brownian motion of such 
a chain in a homogenous flow v(r, t) = ■jit) ■ r may 
be described in a free-draining approximation jj| by a 
Langevin equation 
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ian force with correlations (r](s, t)rj(s' , t')) — 2kT(18(t — 
t')5(s — s'). This equation can be made dimensionlcss 
in terms of reduced variables t = tkT/((L 3 ), s = s/L, 
f = r/L, and L p — L p /L. The linear viscoelasticity 
of a solution of worm-like chains may be characterized 
by either the shear relaxation modulus G(t), which de- 
scribes the stress tr(t) = G(t)[~f + j T ] at time t after 
an infinitesimal step strain 7, or, equivalcntly, by the 
complex modulus G*(u) = iu J^°dt Git)e~ %UJt : , which de- 
scribes the response to a small oscillatory strain. The 
polymer contribution to the moduli per chain, in a di- 
lute solution of c chains per unit volume in a solvent of 
viscosity 77,, is given by a corresponding intrinsic moduli 
\G{t)\ = [G{t)-ri s 5(t)]/c&nd [G*(w)] = [G*{u)-iwri.]/c. 
For worm- like chains, [G{t)\ must have the form [G(t)\ = 
kTG(i, L p ). 

Prior work has identified some relevant time scales and 
provided predictions for G(t) in several limits: 

Rod-like chains (L <c L p ) should behave like rigid rods 
at t > t±, where t± = (L 4 / (kTL p ) is roughly the relax- 
ation time of the longest wavelength bending mode. The 
predicted modulus for dilute rigid rods is 
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Here T is a tension that acts to impose the constraint 
\dr/ds\ = 1, C is a friction coefficient, and rj is a Brown- 



where T roc | = C_L 3 /(72fcT) is a rotational diffusion time. 
The exponential contribution to [G(t)] is due to an en- 
tropic orientational stress caused by an anisotropic dis- 
tribution of rod orientations; it decays by rotational dif- 
fusion. The delta-function contribution arises from the 
longitudinal tension induced in the rods during the step 
deformation; it decays instantaneously after the defor- 
mation. 

In Refs. the authors considered how this behavior 
is modified by the longitudinal compliance of a semiflex- 
ible chain. They calculated the magnitude of changes in 
the end-to-end length of a worm-like chain due to changes 
in the magnitude of transverse fluctuations when the 
chain is subjected to an oscillatory tension at frequency 
ui, and showed that ratio of tension to strain is given by 
a frequency-dependent effective longitudinal modulus 
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at all u 3> rj , where t p = (L 3 /kT. They also pre- 
dicted a macroscopic viscoelastic modulus [G*(u>)] = 
LB(u>)/15 oc (iu:) 3 / 4 jl|||] at very high frequencies, or, 
equivalently, [G(t)] oc t~ 3 / 4 at early times, by assum- 
ing that that the frictional coupling between the chain 
and the solvent must become strong enough at very 
high frequencies to produce an affine longitudinal strain. 
In Refs. the authors considered the dynamics of 

longitudinal relaxation. They showed that longitudinal 
strain propagates along a chain by an anomolous dif- 
fusion with a frequency-dependent diffusivity D(u>) — 
B(u>)/£\\, in which the strain diffuses a distance oc 
y/D(u> = l/t)t oc i 1 / 8 in time t. Both the assumption 
of affine deformation and the predicted t~ 3 ^ 4 decay of 
G(i) must thus fail beyond the time t\\ = (L s /(kTL 5 ) 
required for strain to diffuse the chain length L, and so 
allow signficant longitudinal relaxation. 

This prior work does not predict the behavior of G{t) 
for rod-like chains over a wide range of at intermediate 
times T|| < t < tj_, where relaxation of chain length 
and transverse fluctuations must be coupled. This in- 
terval must rapidly broaden as L <C L p because t\\/tx oc 
(L/L p ) 4 . For L ~ L p , the gaps between tii, tj_, and 
Trod vanish, and so the intermediate regime must dis- 
appear. Coil-like chains [L 3> L p ) are expected || to 
crossover smoothly from G(t) oc t~ 3 l 4 to Rouse-like be- 
havior G(t) oc t^ 1 / 2 at t <~ t p , which is roughly the re- 
laxation time of a bending mode of wavelength L p . 

Our simulations use discrete worm-like chains of N 
beads at positions Ri, . . . , Rjv connected by N — 1 rods 
of fixed length a, with unit tangents = (R»+i — Hi) /a, 

and a bending energy U = — (ft/a) Y^%=2 u » ' u «-i- We 
use a midstep algorithm to compute bead trajectories 
generated by the equation of motion 

fdR 
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Here, Cfc — C a is a bead friction coefficient, F; end = 
—dU/dHi is a bending force, Fj: arlcl is a Langevin noise, 
and F* ens = %Ui — is a constraint force, where 

% is the tension in rod i. The tensions are computed 
by solving H ijTj = • (F i+1 - Fj), where F l = 

Fj end + F™ et + F' and + C&7 • Ri, and H l3 is a tridiagonal 
matrix with Ha = 2 and Hij = — Uj • Uj for i = j ± 1. 
F met = feTgg- In v /det(i?) is a "metric" force that must 
be included in simulations with constrained rod lengths 
to obtain a Boltzmann distribution e -U(ni,...,u N -i)/kT Q f 
rod orientations in thermal equilibrium J^|. 

The modulus [G(i)] is obtained from equilibrium sim- 
ulations (7 = 0) by evaluating the Green-Kubo relation 
[£(*)] = {a xy (t)a xy (0))/kT, where <x = -^VR^ is 
the single-chain stress tensor. The Brownian contribu- 
tion to cr(t) is computed by the method of Refs. MM. 



A wide range of timescales is explored at each value 
of L/L p = Na/L p by using coarser- and finer-grained 
chains to resolve longer and shorter times, respectively. 
Results for chains with the same L/L p but different N 
are collapsed onto master curves of [G(t)] versus i/r ro d, 
where r ro d = C,ba 2 N 3 /72kT . Because initial chain confor- 
mations are chosen from a Boltzmann distribution, be- 
havior at short times can be obtained from short simula- 
tions of fine-grained chains . 

To elucidate the physical origins of stress, we decom- 
pose <r as a sum or = <x ornt + er curv + <r tens — kTI of the 
orientation, curvature, and tension stresses S, where 

Cornt = §fcT(UiUi + UjV-iUAT-x - |l) 

N N-l 
CT curv =-Y / ^F b t Snd +3kTY / (^-^)-^ornt ■ (5) 

i=l z=l 

and <r tens = er cr ornt — cr curv + kTI. We also decompose 
[G(t)] as [G(t)} = [G ornt (t)] + [G curv (t)] + [G te ns(t)], where 
[G a (t)} = (a a , xy (t)a X y(0))/kT, with a = "ornt", "curv", 
or "tens", describes the decay of the stress (<r a {t)) af- 
ter a hypothetical step deformation. er curv arises from 
disturbances of the equilibrium distribution of bending 
mode fluctuations, and was predicted to vanish for rod- 
like chains at times t > t±; er ornt is an analog of the 
orientational stress of a solution of rigid rods; and <r tens 
is the stress arising from longitudinal tension S . 

Fig. (0) shows master curves of [G te ns(i)], [G curv (£)], 
and [Gomt(i)] for chains of reduced length L/L p = 1/8, 
1/4, 1/2. The regions of overlap of results obtained with 
different values of N (which are shown by different col- 
ors), reflect the behavior of a continuous chain, while 
the saturation of [G te ns(t)] an d [G curv (i)] to TV-dependent 
limiting values at small t is due to the discreteness of the 
chains. At long times, t > tj_, the largest contribution 
to [G(t)} is [G or nt(i)], which approaches the exponential 
relaxation predicted for a rigid rod solution. At t ~ t±, 
all three contributions to [G(t)] are comparable. At ear- 
lier times, [G(t)] is dominated by [G te ns(i)]- For the most 
flexible chains shown (L = L p /2), [G te r\s(t)] closely ap- 
proaches the predicted t~ 3 ' 4 asymptote at small t. For 
the two stiffer systems, [Gtens(i)] does not reach this 
asymptote within the accessible range of i, and decays 
more rapidly than t~ 3 l 4 in this range. 

These results are consistent with the prediction of a 
t -3 / 4 decay of [G te ns(^)] below a reduced time T||/r ro d oc 
(L/Lp) 5 that drops rapidly as L/L p decreases, and sug- 
gest the possible existence of a second power law in the 
intermediate time regime tii <C t <C tj_ for L <C L p . 
By postulating the existence of an intermediate power 
law that meets the predicted t~ 3 / 4 asymptote at t ~ tii , 
and that falls to [G ten s(i)] ~ kT at t ~ t±, we obtain 
[Gtens(t)] ~ fcT(t/T ± )- 5 / 4 . The exponent -5/4 agrees 
with the observed slope of log[G ten s(i)] vs. log(f) for the 
stiffest systems shown (L = L p /8), which displays the 
widest intermediate time regime. 
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FIG. 1. Simulation results for [Gtens(i)] (top curve in 
each plot, x), [G CU rv(£)] (middle curve, +), and [G orn t{t)] (bot- 
tom curve, o) vs. t/r rod , with r roc | = ( b N 3 a 2 /(72fcT), for: 
L/Lp = 1/8 and N = 8,16,22,32,46,64,90,128 (top plot) 
L/Lp = 1/4 and N = 8,16,22,32,46,64,90,128 (middle 
plot), and L/L p = 1/2 and N = 16,32,64,90,128,180,256 
(bottom plot). In each plot, [G or nt(i)] is shown only for 
a few small values of N. The long-dashed red line is the 
prediction [G ro d(t)} = §fcTe -t/Trod for rigid rods at t > 0. 
Short-dashed red lines with slopes of —3/4 and —5/4 are 
asymptotes Eqs. ([li]) and (|l3|). The solid red line is the pre- 
dicted [Gtens(i)] obtained by Fourier transforming Eq. (Jl l|) . 



We now present a theory of the longitudinal dynamics 
of a rod-like chain that predicts the observed t~ 5 / 4 de- 
cay of [Gtens(i)] at intermediate times for L -C L p . Our 
analysis resembles one given previously to describe lon- 
gitudinal relaxation of entangled chains S . For rod- like 
chains, we may expand r(s, t) around a rod-like reference 
state as r(s, t) — r\\ (s, t)n(t) + h(s, t), where h(s, t) satis- 
fies h(s, i)-n(t) = 0, and n(i) is a unit vector that rotates 
with the flow like a non-Brownian rigid rod: ri = P • 7 • n, 
where P = I — nn. Linearizing Eq. ([!]) then yields longi- 
tudinal and transverse equations 



C ^~ r ll7:nn} = — + m 
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These equations are coupled by the tension T, which is 

chosen to satisfy the constraint |§§| 2 = l^fP + lfjl 2 = 1- 
It is convenient to introduce a longitudinal strain field 



e(s) 



dr it (s) 



(-rp-)eq, where (■ • -) eq denotes a thermal 



ds \ ds /eq' vvii ^ \ /eq 

equilibrium average. Combining this definition with the 
constraint and expanding to lowest order in |<9h/9s| 2 
yields an approximate expression of e(s) in terms of 
h(s), e(s) ~ -l{|^i|2_ ( |^|2 >eq }. This expres _ 

sion for e(s) and Eq. (Q) were used in Refs dH to cal- 
culate the linear response of the spatial average strain 
(e(u))) = J ds(e(s, cj)) /L to a spatially uniform oscillating 
tension T(uS) at frequency u> (where functions of u> de- 
note Fourier amplitudes), yielding an effective extension 
modulus B(uj) = T(w)/(e(u>)) that is given by Eq. (||) at 
w > tI 1 ! li>i and b y a st atic value B(0) ~ fcTi^/L 3 at 
- • 

A modified diffusion equation for the strain field may 
be obtained by taking the thermal average of Eq. (^) , dif- 
ferentiating with respect to s, Fourier transforming with 
respect to t, and setting (T(s,uj)} — B(u>){e(s, uSj). This 
yields 



(e(s,u))) ~ 1^7(0;) 



nn 



(8) 



where B(cu)/( is an effective diffusivity and 7(w) is the 
amplitude of an oscillatory strain tensor. Hereafter n(i) 
is approximated by its time average over one period of 
oscillation. Eq. (||), with (e(s,w)) = at the chain ends, 
has the solution 



1 - 



cosh(A(w)(2s - 1)) 



cosh(A(u/)) 



7(0;) : nn 



(9) 



where = (iu(L 2 /4B(w))Va = ( iW r [ |/2 11 ) 1 /8. 

The tension stress of rod-like chains subjected to an 
infinitesimal oscillatory strain 7(0;) is given by 



<T te 



5 M 



/ ds (T(s,u) nn) 
Jo 



(10) 
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where (■ • ■) denotes an average over both weak fluctua- 
tions and overall rod orientations. Combining Eq. JI(i| ) 
with Eq. (H) for the strain along a rod of orientation 
n and averaging over random rod orientations yields a 



stress <T t ens(w) = [G*, 

lus 



.(cj)][7((jj) +7 (&>)], with a modu- 



1 - 



tanh(A(cj)) 



A(w) 



(11) 



This prediction has the following limiting behaviors: At 
frequencies lu ^> TiT" 1 , Eq. ( |ll| ) reduces to the high- 
frequency asymptote [G* en5 (o;)] ~ LB(u>)/15 found pre- 
viously (iU|] . Fourier transforming this asymptote yields 
a relaxation modulus 



lim [G tens (i)] 

*<T|| 



G r 



-3/4 



where d = 2 3 / 4 /[15r(i)] = 0.0309. At intermediate 
frequencies r. 



> to > t7 



where A(w) <C 1, ex- 
panding Eq. (11) in powers of A(u>) yields a modulus 



180 



fcT 
1800 2 3 / 4 



(iWT ± ) 5 / 4 



which in- 



cludes a dominant contribution of order iu>, whose prefac- 
tor is identical to that found for rigid rods, and a first cor- 
rection proportional to (icu) 5 ^ 4 . This yields a loss mod- 
ulus [G"(o;)] oc u> (like rigid rods) but a storage modu- 
lus [G'(w)] oc lu 5 / 4 (unlike rigid rods) at these frequen- 
cies. Upon transforming this intermediate asymptote, 
the term proportional to iuj yields an apparent (5-function 
contribution to [G(i)] (as for rigid rods), and so [G ten s(*)] 
is instead dominated at n <C t <C tj_ by the transform 
of the term of proportional to (iuj) 5 ^ 4 , which yields 



lim [G ten5 

rii<t<TL L 



(*)] 



< / I — 



-5/4 
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where G 2 = l/[2 3 / 4 7200r(f )] = 0.0000674. At w < tJ 1 
or t > t_l, Eq. (j|) for B{cu) becomes inapplicable, but 
[Gtens(i)] & l so becomes small compared to [G wt(t)]- 

The predictions of [Gtens(i)] shown in Fig. (|l|) were 
obtained by Fourier transforming Eq. ( |ll"l) numerically. 
They agree with the simulation results for [G te ns(i)] at all 
t J$ tj_, and accurately describe not just the power law 
regimes, but the broad crossovers between them. Re- 
markably, the theory remains accurate for L = L p /2, 
despite the assumption of a nearly straight chain. 

Our derivation of Eq. (||) explicitly assumes a local 
proportionality of the tension and strain at each point 
on the chain, with (e(s,u>)) — B^ 1 (uj)(T(s,uj)), rather 
than allowing for a spatially nonlocal response of the form 
(e(s,uj)} = Jds' B^ 1 (s,s',uj)(T(s',uj)). To examine this 
approximation, we calculated the nonlocal compliance 
B^ 1 (s, s' , uj). We find that the range of nonlocality is of 
the order of the wavelength £j_(w) = (cj^/kTLp)~ 1,>4 of 
the bending mode with frequency u>, and that the strain 



predicted by Eq. M) varies slowly over lengths of order 



£_i_(w) for all w 



> 



This justifies our local compliance 



approximation for all u> > r 7 . 

A conceptually simple, analytically solvable, and ac- 
curate model of the dominant contribution to [G(t)] at 
times t < t±, valid for all L < L p , is thus obtained 
by treating the inextensible worm-like chain as an effec- 
tively extensible rod with a frequency-dependent longi- 
tudinal modulus given by Eq. @. At later times, [G(t)] 
is dominated by [G om t(^)]i which mimics the behavior of 
a solution of rods. The simulations show that the curva- 
ture stress never dominates [G(t)] in such solutions. A 
useful global approximation for [G(i)] for rod-like chains 
may thus be obtained simply by replacing the 5-function 
in Eq. (||) by our result for [G ten sW]- O ur results confirm 
that [G(t)} initially decays as t~ 3 / 4 , but also show that, 
(12) when L < L p , this behavior is observable only below a 
time proportional to rn that drops rapidly with decreas- 
ing L/Lp to values inaccessible to either simulation or 
experiment. Therefore, measurements of the viscoelastic 
modulus of dilute solutions of rod-like chains at practi- 
cally attainable high-frequencies may often probe either 
the t~ 5 / 4 regime identified here, instead of the initial 
£—3/4 re gj mej or the broad — but calculable — crossover 
between them. 
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